Optimisation en architecture

Kayané Robach

Magistère d'Economiste Statisticien

11 juin 2021

Le but de ce projet est de construire et résoudre un problème d’optimisation appliqué à l’architecture. Initialement l'objectif était de proposer une alternative aux logements sociaux disposés en blocs. Pour cela nous définirons des contraintes à respecter, à savoir, minimiser le vis à vis, minimiser les surfaces de contact entre habitations, maximiser le champ visuel.

Nous nous intéressons donc au problème suivant :

Comment construire une structure agréable à vivre en utilisant un recuit simulé ?

Nous verrons dans un premier temps comment répondre à ce problème en 2D en utilisant la méthode du recuit simulé. Puis nous appliquerons ce processus en 3D pour répondre à notre problématique. Le choix de cette méthode mathématique s'est imposé au vue du cadre discret du problème et de la complexité de l'analyse. En effet je n'ai trouvé aucune autre méthode pertinente pour répondre aux objectifs choisis.

Algorithme du recuit simulé

Soit $E$, l'espace d'état, fini et $H : E \longrightarrow \mathbb{R}$, la fonction objectif, une fonction dont on cherche un minimum.

On considère une mesure de Gibbs paramétrée par $T > 0$ de la forme $$\mu_T(x) = \frac{1}{Z_T} exp^{- \frac{H(x)}{T}}$$ où $Z_T$ est la constante de normalisation souvent appelée fonction de partition. Lorsque $T \longrightarrow 0$ la mesure $\mu_T$ se concentre sur les cubes de E proche des minimums de H ; on considère donc une suite $(T_n)_n$ telle que $\underset{n \longrightarrow \infty}{lim} T_n = 0$ et l'algorithme de Hastings-Metropolis afin trouver la solution minimum à notre problème.

L'algorithme de Hastings-Metropolis permet de construire une chaîne de Markov de probabilité invariante $\mu$ donnée ; en principe cet algorithme est utilisé pour simuler une loi (grâce à une chaine de Markov) qui est difficile à simuler.

L'algorithme construit de manière itérative une suite $(x_n)_n$ d'éléments de E. Le terme $x_n$ de la suite $(x_n)_n$ désigne le cube du chemin $(x_n)_n$ envisagé que nous sommes en train de construire et qui parcoure E dans un but précis (dépendant de chaque problème). Pour chaque nouvelle itération cherchons un cube suivant qui nous rapprochera de la solution du problème dont il est question. Nous sélectionnons donc un nouveau chemin candidat $(y_n)_n$ puis nous le testons afin de décider si nous le choisissons ou non comme étape suivante de notre chemin solution.

Une matrice $P$, dite de sélection, est définie afin de parcourir l'espace d'état E. Notre but est de construire une chaîne de Markov $(X_n)_n$ de probabilité invariante $\mu$. La matrice de sélection permet de sélectionner des candidats $y = (y_n)_n$ voisins de notre chemin actuel $x = (x_n)_n$ pour construire le chemin solution représenté par $(X_n)_n$. Ainsi $\forall x, y \in E, P(x,y) = \mathbb{P}(Y_n^{(x)} = y)$ défini la probabilité que l'on choisisse le chemin y en tant que candidat (la variable aléatoire du candidat étant $(Y_n)_n$) à partir du chemin actuel $x$.

Une fonction $\rho$ est ensuite définie à partir de $P$ et de $\mu$ par : $$\forall x,y \in E, \rho(x,y) = \frac{\mu(y) P(y,x)}{\mu(x) P(x,y)} \wedge 1$$ $\rho$ c'est la probabilité que l'on accepte le candidat $y$ selectionné par $P$.

A partir de toutes les définitions présentées nous pouvons ensuite construire la matrice stochastique Q comme la probabilité de procéder à la transition de x à y à l'étape n de la construction de notre suite solution. En termes mathématiques : $$Q(x,y) = \mathbb{P}[X_{n+1} = y | X_n = x] = \left\{ \begin{array}{ll} P(x,y) \rho(x,y) & \mbox{si } x \neq y \\ 1 - \sum_{z \neq x} Q(x,z) & \mbox{sinon.} \end{array} \right.$$

Sous les hypothèses adhoc, la chaîne de Markov $(X_n)_n$ (de transition Q) converge en loi vers $\mu$.

Pour résumer, on choisit un voisin $y$ de $x$ avec probabilité $P(x,y)$ (donnée par la matrice de sélection) puis, avec probabilité $\rho(x,y)$ on accepte ce voisin, auquel cas on pose $x_{n+1} = y$ ou on le refuse et on pose $x_{n+1} = x$

Concrètement la définition de notre chaîne de Markov est la suivante : (nous avons besoin de $(U_n)_n \approx \mathcal{U}_{[0,1]}$ pour construire $(X_n)_n$ rigoureusement) $$X_{n+1} = Y_{n+1}^{X_n} \mathbb{1}_{\{U_{n+1} \leq \rho(X_n,Y_{n+1}^{X_n})\}} + X_n \mathbb{1}_{\{U_{n+1} > \rho(X_n,Y_{n+1}^{X_n})\}}$$

Dans le cadre du recuit simulé la probabilité de procéder à la transition de x à y à l'étape n de la construction de notre suite solution : $Q$ est la suivante : $$Q_n(x,y) = \mathbb{P}[X_{n+1} = y | X_n = x] = \left\{ \begin{array}{ll} P(x,y) exp^{- \frac{(H(y) - H(x))_{+}}{T_n}} & \mbox{si } x \neq y \\ 1 - \sum_{z \neq x} Q_n(x,z) & \mbox{sinon.} \end{array} \right.$$

Voici un résumé schématique du fonctionnement du recuit simulé :

Problème du petit lotissement

Sur un terrain carré nous voulons construire des maisons vitrées. Une contrainte doit être vérifiée : il doit y avoir le moins de vis à vis possible. Pour cela la première étape est d'éloigner les maisons au maximum les unes des autres.

Pour répondre à ce problème d'optimisation nous utiliserons l'algorithme du recuit simulé expliqué ci dessous.

Nous travaillons pour le moment en dimension 2.

Modélisation de notre problème du petit lotissement

L'espace d'état E ici est une grille que l'on représentera par une matrice $k$ x $k$. Nous considèrerons $p$ cubes de cette grille défini chacun par 1 entier difféfent de 0 comme étant les maisons. La matrice "terrain" contient aisni des 0 pour les espaces libres et des entiers allant de 1 à $p$ pour les habitations.

Soit p maisons représentées par un vecteur $M = (M_1, M_2, M_3, ..., M_p)$ ; chaque entier représentant une maison se voit attribuer des coordonnées dans la matrice en 2 dimensions $\in [0,1,2, ..., k] \times [0,1,2, ..., k]$. Nous notons $d$ la distance euclidienne entre 2 maisons $M_i$ et $M_j$ $i,j \in {1,2,3, ..., p}$ : $d(M_i, M_j)$.

Notre algorithme devra retourner un grille remplie avec des maisons (donc une matrice de 0 avec un 1, un 2, un 3, ... et un $p$ pour représenter chaque maison dans la grille).

On souhaite maximiser la distance entre les maisons dans notre petit lotissement. On a donc un vecteur de distances $( d(M_1, M_2), d(M_1, M_3), d(M_2, M_3), ... )$ à maximiser.

H, la fonction objectif que nous allons minimiser, peut être donc être définie dans un premier temps comme l'opposé de la norme du vecteur des distances. Avec $p$ = 3 on a : $$H(M) = - ||(d(M_1, M_2), d(M_1, M_3), d(M_2, M_3))||$$

Code de la distance

On utilisera la distance euclidienne.

Code de H

Pour commencer, la fonction objectif est l'opposé de la norme du vecteur des distances. Si on veut maximiser les distances, on veut donc maximiser la norme du vecteur des distances. Cela revient à minimiser l'opposé.

Recuit simulé

Concrètement, nous avons une matrice avec des emplacements occupés et des emplacements libres. Nous choisissons une maison à déplacer. Nous cherchons les emplacements voisins de cette maison qui sont libres et nous sélectionnons un de ces espaces non occupés.

$ \left\{ \begin{array}{ll} \mbox{si } H \searrow & \mbox{nous déplaçons la maison }\\ \mbox{sinon } & \mbox{nous déplaçons la maison de temps en temps } \end{array} \right. $

Nous souhaitons explorer de temps en temps les alentours même si H ne diminue pas (d'où l'idée d'algorithme stochastique). Cela nous permettrait de sortir d'un minima local si nous nous retrouvions coincés.

Nous reviendrons plus en détails sur ce processus plus tard.

Définissons la matrice de sélection $P$ qui nous permettra de parcourir l'expace ainsi qu'une suite $(T_n)_n$ telle que $\underset{n \longrightarrow \infty}{lim} T_n = 0$.

Pour cela nous allons donner une matrice initiale à l’algorithme (représentant ici notre terrain avec quelques maisons). Nous choisissons une maison aléatoirement, puis avec une certaine probabilité nous déplaçons cette maison.

Code de la suite T

Un théoreme énonce que si la matrice de sélection P est symétrique irréductible et si cette matrice vérifie la condition de Doeblin (faisant intervenir les chaînes de Markov et la théorie des probabilités) et si $$\forall n \geq 1, T_n = \frac{\gamma}{log(n)}$$ avec $\gamma > \underset{x, y}{max}(H(y) - H(x)) \mathbb{1}_{\{P(x,y)>0\}}$

Alors pour toute loi initiale $\nu_0$ de $X_0$ on a : $$\underset{n \longrightarrow \infty}{lim} ||\nu_n - \mu_n|| = 0$$ donc en particulier, $\underset{n \longrightarrow \infty}{lim} P( H(X_n)>argmin H ) = 0$

Code de la matrice de sélection

Pour le moment, étant donné une combinaison de position de maisons (un trajet dans le cas du problème du voyageur), on tire une maison aléatoirement, et on bouge cette maison dans un des espaces libres autour d'elle avec une certaine probabilité. Autour d'elle signifie dans un premier temps : à une norme 1 vallant 1.

Pour construire $P$ comme on le souhaite il faut d'abord une fonction qui identifie les espaces libres autour d'une maison.

Pour cela, nous créons une matrice d'indices (ou plutot 2 pour les indices lignes et les indices colonnes) que nous centrons avec les coordonnées de notre point.

Prenons un exemple : nous regardons le point qui devrait se trouver dans le carré rouge (point d'indice (2,2) car on commence à 0).

Nous construisons une matrice d'indices lignes et une matrice d'indices colonnes.

Nous centrons les matrices en soustrayant l'indice ligne à la matrice d'indices lignes et l'indice colonne à la matrice d'indices colonnes.

Nous calculons ensuite la somme des valeurs absolues des coordonnées des 2 matrices d'indices. Dans la matrice obtenue, les éléments $\leq$ 1 sont les voisins de notre point.

Nous pouvons ensuite construire $P$. Après sélection d'une maison et identification de l'espace libre autour de cette dernière, nous choisissons aléatoirement un emplacement parmi les emplacements libres et nous déplaçons la maison.

Revenons à notre problème global, nous cherchons à minimiser une fonction H que nous connaissons car nous l'avons défini mais dont nous ne connaissons ni la forme ni le fond. Pour éviter les minima locaux nous faisons donc fluctuer notre avancée de manière aléatoire (en déplacant aléatoirement un cube à une place libre choisie aussi aléatoirement). C'est ainsi que nous seront capable d'aboutir au minimum de H (nous l'espérons).

Rappelons la définition de la marice de transition de notre problème de recuit simulé.

$$Q(x,y) = \mathbb{P}[X_{n+1} = y | X_n = x] = \left\{ \begin{array}{ll} P(x,y) \rho(x,y) & \mbox{si } x \neq y \\ 1 - \sum_{z \neq x} Q(x,z) & \mbox{sinon.} \end{array} \right.$$$$\rho(x,y) = \frac{\mu(y) P(y,x)}{\mu(x) P(x,y)} \wedge 1$$$$Q(x,y) = \mathbb{P}[X_{n+1} = y | X_n = x] = \left\{ \begin{array}{ll} P(x,y) (\frac{\mu(y) P(y,x)}{\mu(x) P(x,y)} \wedge 1) & \mbox{si } x \neq y \\ 1 - \sum_{z \neq x} Q(x,z) & \mbox{sinon.} \end{array} \right.$$$$\mu_{T_n}(x) = \frac{1}{Z_{T_n}} exp^{- \frac{H(x)}{T_n}}$$$$Q(x,y) = \mathbb{P}[X_{n+1} = y | X_n = x] = \left\{ \begin{array}{ll} P(x,y) (\frac{\frac{1}{Z_{T_n}} exp^{- \frac{H(y)}{T_n}} P(y,x)}{\frac{1}{Z_{T_n}} exp^{- \frac{H(x)}{T_n}} P(x,y)} \wedge 1) & \mbox{si } x \neq y \\ 1 - \sum_{z \neq x} Q(x,z) & \mbox{sinon.} \end{array} \right.$$

$P$ est supposée symétrique.

$$Q(x,y) = \mathbb{P}[X_{n+1} = y | X_n = x] = \left\{ \begin{array}{ll} P(x,y) ( \frac{ exp^{ - \frac{H(y)}{T_n}} }{ exp^{- \frac{H(x)}{T_n}} } \wedge 1) & \mbox{si } x \neq y \\ 1 - \sum_{z \neq x} Q(x,z) & \mbox{sinon.} \end{array} \right.$$$$Q(x,y) = \mathbb{P}[X_{n+1} = y | X_n = x] = \left\{ \begin{array}{ll} P(x,y) (exp^{- \frac{H(y)}{T_n} + \frac{H(x)}{T_n}} \wedge 1) & \mbox{si } x \neq y \\ 1 - \sum_{z \neq x} Q(x,z) & \mbox{sinon.} \end{array} \right.$$$$ exp^{- \frac{H(y)}{T_n} + \frac{H(x)}{T_n} } < 1 \iff exp^{- \frac{H(y)}{T_n} + \frac{H(x)}{T_n} } < exp^{0} \iff \frac{H(x)}{T_n} - \frac{H(y)}{T_n} < 0 $$$$ \implies \left\{ \begin{array}{ll} exp^{- \frac{H(y)}{T_n} + \frac{H(x)}{T_n}} \wedge 1 = 1 & \mbox{si } \frac{H(x)}{T_n} - \frac{H(y)}{T_n} \geq 0 \\ exp^{- \frac{H(y)}{T_n} + \frac{H(x)}{T_n}} \wedge 1 = exp^{- \frac{H(y)}{T_n} + \frac{H(x)}{T_n}} < 1 & \mbox{sinon.} \end{array} \right. $$$$ \iff \left\{ \begin{array}{ll} exp^{- \frac{H(y)}{T_n} + \frac{H(x)}{T_n}} \wedge 1 = exp^{0} = 1 & \mbox{si } (H(y) - H(x))_{+} = 0 \geq 0 \\ exp^{- \frac{H(y)}{T_n} + \frac{H(x)}{T_n}} \wedge 1 = exp^{- \frac{H(y)}{T_n} + \frac{H(x)}{T_n}} < 1 & \mbox{sinon.} \end{array} \right. $$$$Q_n(x,y) = \mathbb{P}[X_{n+1} = y | X_n = x] = \left\{ \begin{array}{ll} P(x,y) exp^{- \frac{(H(y) - H(x))_{+}}{T_n}} & \mbox{si } x \neq y \\ 1 - \sum_{z \neq x} Q_n(x,z) & \mbox{sinon.} \end{array} \right.$$

Concrètement que doit-il se passer ?

On a notre matrice avec une certaine combinaison d'emplacements occupés et d'emplacements libres. On choisi aléatoirement une maison à déplacer (grâce à la matrice de sélection $P$) et, on choisi le nouvel emplacement qu'elle va éventuellement occuper (parmi les emplacements voisins de cette maison). La norme utilisée a ici un rôle très important dans la définition des espaces voisins disponibles et/ou occupés. Avec une certaine probabilité (définie par $\rho$) nous déplaçons cette maison. Comment ?

Nous avons là défini l'éventuel nouvel emplacement de notre maison. Nous comparons alors la valeur de notre fonction objectif $H$ évaluée à l'ancienne matrice puis à la nouvelle (après éventuel délacement d'une des maisons).

Si $H(new\_matrix) < H(matrix)$, alors nous déplaçons la maison.

Si $\mathcal{U}_{[0,1]}$ < $exp^{ \frac{H(matrix) - H(new\_matrix)}{T(n)} }$, alors nous déplaçons la maison.

L'Algorithme de Metropolis

Let's go !

C'est parti, nous lançons l'algorithme sur une matrice 10x10 en 2D.

Un problème apparaît. Notre fonction H étant définie comme elle l'est, c'est à dire, en prenant en compte la distance moyenne des maisons aux autres, l'algorithme n'a pas d'intérêt à changer une combinaison comme celle qui apparaît à la fin. En effet, si on déplace une des maisons qui est acollée à une autre, notre fonction objectif n'est pas modifiée.

Nous allons donc changer notre fonction objectif.

Petit retour en arrière pour redéfinir H.

Code de H

La fonction objectif est maintenant la somme des inverses des distances de chaque maison aux autres. Si on veut maximiser les distances, on veut donc minimiser les inverses des distances. Cela permet de mettre plus de poids sur les distances courtes et ainsi de faire en sorte que chaque maison s'isole des autres.

Voyons si nous avons réussi améliorer notre algorithme :

Regardons ce que ca donne quand la moitié de la surface est occupée :

Problème du HAV, l'Habitation Agréable à Vivre

Nous avons vu le problème du petit lotissement ; fini les plaisanteries, nous passons désormais en 3D. Sur un terrain carré (ou plutôt cubique car nous pouvons considérer la hauteur), nous voulons construire un batiment.

Ce sera une structure où chacun des cubes qui la composent représente une habitation. Comme tout à l'heure nous imaginons des petites habitations vitrées. Mais attention, nous changeons les contraintes : il doit toujours y avoir le moins de vis à vis possible et nous voulons que les habitations ne forment pas un block, nous souhatons une structure plus épurée. Nous allons définir proprement toutes ces contraintes par la suite.

Nous travaillons désormais en 3D ; dans notre matrice initiale, tous les cubes sont au sols.

Voici notre matrice initiale

Pour répondre à nos besoins nous définissons quelques fonctions utiles.

Nous avons besoin de 4 fonctions pour cela ; get_matrix_indexes fait, tout comme en 2D, une matrice d'indices, que cette fois nous représentons en tuple (x, y, z). centered_matrix_norm existe en 2 exemplaire, un pour la norm 1 et un pour la norm 2. Tout comme en 2D encore, ces fonctions centrent la matrice des indices à partir des coordonnées d'un cube spécifique afin d'identifier les voisins du cube dont il est question. Enfin get_location_around sert à la fois à identifier les emplacements libres et les emplacements occupées autour du cube passé en paramètre.

Grâce aux fonctions ci dessus nous allons définir nos objectifs. Nous voulons désormais minimiser les surfaces en contact et maximiser le champ visuel.

Premièrement nous voulons minimiser les surfaces en contact entre cubes, donc H prend en compte la moyenne des surface en contact par cube. Ensuite nous voulons maximiser le champ visuel, donc minimiser l'inverse du champ visuel. Rappelez vous comment nous avons défini le champ visuel : il est le total des cases séparant la face d'un cube et le prochain cube rencontré sur le même axe (pour une même hauteur). Pour être plus explicit, imagnez un cube-appartement traversant (avec des fenêtres sur 2 faces opposées), vous vous tenez au milieu de la pièce et regardez une fenêtre puis vous comptez l'espace entre votre appartement-cube et le suivant qui se trouve dans votre champ de vision ; en sommant ces valeurs pour chacune des fenêtres nous obtenons ce que nous avons appelé ici le champ visuel.

Mais il y a 2 champs visuels, le Nord-Sud et le Est-Ouest. A chaque fois, nous avons donc décidé de considérer comme le champ visuel à maximiser, celui qui est le plus grand.

Ainsi H est la somme de :

Code de H

C'est parti pour définir la matrice de passage P cette fois en 3D.

Comment faire ?

Nous piochons un cube et nous identifions les espaces libres et occupés qu'il a autour de lui. (Et s'il n'y en a pas nous en piochons un autre).

Là, nous allons proposer une nouvelle place pour le cube sélectionné, et nous allons calculer le nombre de groupes de cubes qui seraient alors identifiés (grâce à la matrice d'affinité), si nous déplacions le cube. Rappelez vous, nous souhaitons construire une structure donc nous voulons identifier 1 seul groupe car il faut que nos cubes soient tous "liés" entre eux.

Donc, tant que la place proposée pour le cube sépare notre structure nous en choisissons une autre.

Plusieurs paramètres sont définis pour controller notre algorithme :

Lançons l'algorithme :

Nous afficherons les étapes de la "construction" lorsque l'algorithme approche une impasse.

Avez vous identifié le problème ?

Il arrive qu'un cube que l'on déplace sépare la structure, la faisant passer d'un seul et même groupe à 3 groupes distincts de cubes. Le cube que l'on déplace à ce moment ne peut jamais palier à ce problème sans revenir à sa position initiale (d'où on le déplace justement).

Nous avons donc plusieurs possibilité devant nous :

Jusqu'à présent nous utilisions la norme 1 égale à 1 pour définir une proximité entre cubes.

Testons tout d'abord avec la norme 2 afin de considérer comme "en contact" des cubes qui se trouvent en diagonale.

C'est difficile de garder une structure en 1 morceau avec cette solution...

Maintenant, toujours avec la norme 2, nous commençons avec une plus grosse quantité de cubes. En changeant ainsi la matrice initiale et la quantité de cubes nous espérons aboutir à un résultat plus convaincant.

Visiblement, l'algorithme ne rencontre pas d'erreur avec un plus grand nombre de cubes. Néanmoins, comparé aux situations précédentes nous pouvons voir sur ce graphique d'évolution de la fonction H que le minimum n'est pas atteint. L'algorithme met beaucoup de temps à s'éxécuter, il faudrait augmenter le nombre d'itérations. Le même problème persiste, avec la norme 2 il est impossible d'obtenir une structure viable.

Nous allons donc tenter d’utiliser la norme 2 uniquement pour identifier le nombre de groupes de cubes (des cubes en diagonale sont donc considérés comme étant en lien). Nous gardons la norme 1 pour le reste, c’est à dire la localisation des voisins des cubes, les surfaces en contact, l’identification des places libres autour... La norme 2 interviendra uniquement dans le calcul de la matrice de liens pour définir les affinités entre cubes. (Pour cela nous définissons à nouveau la fonction "create_affinity_matrix").

Encore une fois, le code est très long et rencontre une erreur avant la fin...

Pour conclure, nous avons obtenu la meilleure combinaison avec la norme 1 et un "grand" nombre de cubes. Voici une matrice initiale accompagné de 2 résultats obtenus lors de différentes executions du code.

Ces résultats finaux correspondent bien à nos objectifs mêmes s'ils ne représentent pas le résultat optimal atendu minimisant la fonction objectif.

Finalement, pour répondre à notre problématique initial, il semble important d'optimiser son code car la procédure du recuit simulé peut être très longue. De plus, comme nous l'avons vu, la définition de la fonction objectif à minimiser est une étape importante et complexe de ce projet. Globalement sur les 2 résultats finaux présentés ci-dessus, une structure arborescente apparaît. Un problème éventuel serait le fait que des cubes s'accrochent côte à côte sans qu'aucun cube ne soit présent en dessous ; ce qui n'est pas possible, il faudrait donc prendre en compte cet aspect là dans la construction si nous voulions continuer ce projet.

Néanmoins, les résultats obtenus sont en accord avec les contraintes définies. L'objectif de ce projet est ainsi en parti atteint et laisse place à de nouvelles perspectives pour prolonger et approfondir cette recherche.